function [vals, labels, VarNames] = makeArrayFromFullCartesianCSV(fileName, NumHeaderCols, varargin)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function reads a CSV file where the first N columns are "headers" (user-id, week-id, product-id, ...),
	% and all the other column are numeric values. Each combination of id's should be unique (there should be at most
	% one row for a given combination of id's), and there should be one line for each such combination (i.e. the CSV
	% contains the full cartesian product of id's).
	%
	% The file is set up as output from function fillCartesianProduct in SQLite (defined in SQLite/utilities.sh).
	% For example:
	%
	% ID1,ID2,ID3,val1,val2
	% A,1,a,3,5
	% A,1,b,4,11
	% A,2,a,5,23
	% A,2,b,6,12
	% B,1,a,7,7
	% B,1,b,8,3
	% B,2,a,9,0
	% B,2,b,10,7
	% 
	% In this example: NumHeaderCols should be set to 3:
	% values = makeArrayFromFullCartesianCSV(myfile.csv, 3);
	%
	% or rather (because of non-numeric IDs):
	% values = makeArrayFromFullCartesianCSV(myfile.csv, 3, {{'char', 'double', 'char'}});
	%
	% The function reads the CSV file as a table, discards the first 3 columns, converts the table into an array,
	% and resizes its dimensions as: NumIds1 x NumIds2 x NumIds3 x 2  (because there is val1 and val2)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% fileName:								
	% NumHeaderCols:			integer (number of columns that are "labels" - not "numeric" data)string
	% varargin{1}=headerTypes:	cell(NumHeaderCols, 1)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% values:					NumDistinctHeadersN x ... NumDistinctHeaders2 x NumDistinctHeaders1 x NumValues
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	opts = detectImportOptions(fileName);
	NumCols = length(opts.VariableTypes);
	if nargin >= 3
		headerTypes = varargin{1};
		for jj = 1:NumHeaderCols
			opts.VariableTypes{jj} = headerTypes{jj};
		end
	end
	for jj = NumHeaderCols+1:NumCols
		opts.VariableTypes{jj} = 'double';
	end
	VarNames{1} = opts.VariableNames(1:NumHeaderCols);
	VarNames{2} = opts.VariableNames(NumHeaderCols+1:NumCols);

	M = readtable(fileName, opts);

	NumRows = size(M, 1);
	NumVars = NumCols - NumHeaderCols;
	
	mydims = zeros(1, NumHeaderCols);
	for jj=1:NumHeaderCols
		mydims(jj) = size(unique(M(:,jj)), 1);
	end
	mydims = [fliplr(mydims) NumVars];

	vals = M(:,NumHeaderCols+1:end);
	vals = table2array(vals);
	vals = reshape(vals, mydims);
	
	warning('PLEASE MAKE SURE THAT THE ORDERING IS CONSISTENT!!!');
	
	% Get the labels
	distinctHeaderValues = cell(NumHeaderCols, 1);
	for jj=1:NumHeaderCols
		labels{jj} = sort(table2array(unique(M(:,jj))));
	end
end
